Chemical potential in 3D lattice¶

Contact oleg.malanyuk@epfl.ch if you have any problems or feedbacks to this notebook

Introduction¶

This notebook gives a simple illustration to the variation of chemical potential for Fermions in 3D, following the Fermi-Dirac statistics in finite temperature.

In [1]:
import matplotlib.pyplot as plt
import numpy as np
import math as m
import scipy
from ipywidgets import interactive as jupyter_interactive
from ipywidgets import interact
from ipywidgets import Layout
from scipy import integrate
import ipywidgets as ipw

from bokeh.models import ColumnDataSource
from bokeh.io import push_notebook, output_notebook, show
from bokeh.layouts import row, column
from bokeh.plotting import figure
output_notebook()

%matplotlib inline

# Basic libraries
Loading BokehJS ...

In 3D, we have the density of states $\rho \propto \sqrt{E}$. At zero temperature, the distribution of electrons just follows $$g(E)=a \sqrt{E}$$ and the Fermi energy is determined by the total occupied number $N_{occ}$ :$E_f$ is the energy such that $$N_{occ}=\int_{0}^{E_f} g(E) dE$$

In finite temperature, however, one must also consider the Fermi-dirac distribution of $$f(E)=\frac{1}{e^{(E-\mu)/kbT}+1}$$ when calculating the total occupation number. The expression is then $$ N_{occ}=\int_{0}^{\infty} f(E)g(E) dE$$

The chemical potential $\mu$ must be consistent with the occupation number. As shown in the following, the chemical potential changes at finite temperature.

The above plots give a simple illustration of variation of the chemical potential at finite temperature. At finite temperature, the simple calculation with zero-temperature $E_f$ is no longer consistent with the occupation number. The calculation of the varied fermi energy is shown in chapter 2.2.3 with some simple approximation.

Alternatively, we can also calculate the variation of chemical potential with numerical methods (trying to fit the chemical potential such that the integral at finite T gives same electron density as at T = 0K). In the content of numerical evaluation, the consistency problem of chemical potantial and the occupation number can be solved through iteration.

In [2]:
#Initial Plot, kT taken as 0.1Ef

kt = 0.1
ef = 1
kBT = 0.001 # Boltzmann constant measured interms of chemical potential
E_min = 0.0001
E_max = 3.0*ef
factor = 1.0
N = 100
n0=50
N1=20
# Set up the variable space for 1D plot
En_mu = np.linspace(E_min, E_max, N) #Energy measured in units of chemical potential
En_l=np.linspace(E_min,ef,n0)
En_r=np.linspace(ef,E_max,n0)
ef1=ef
nocc=2.0/3.0*(ef**(3.0/2.0))

def nfg(e0):
    v,err=scipy.integrate.quad(cf,0,e0)
    return v
def fg(e0):
    f1=lambda x: np.sqrt(x)
    v,err=scipy.integrate.quad(f1,0,e0)
    return v

num_s=np.linspace(1,N1,N1)
ef1_s=np.linspace(ef,ef,N1)
for num in range(1,N1):
    cf=lambda x:np.sqrt(x)/(np.exp( (factor*x-ef1)/kt )+1) 
    c_f = cf(En_mu)  #np.sqrt(En_mu)/(np.exp( (factor*En_mu-1)/kt )+1) # Occupation number according to F-D statistics
    c_df = np.gradient(c_f) # Take the first derivative of the F-D distribution.


    #print(nfg(0.500))
    n_e=(1.5*nfg(E_max))**(2.0/3.0)

    ef1_s[num]=ef1
    ef1=ef1-0.2*(n_e-ef)

ff0=lambda x:1.0/(np.exp( (factor*x-ef)/kt )+1)
cf0=lambda x:np.sqrt(x)/(np.exp( (factor*x-ef)/kt )+1) 
ff=lambda x:1.0/(np.exp( (factor*x-ef1)/kt )+1)
cf=lambda x:np.sqrt(x)/(np.exp( (factor*x-ef1)/kt )+1) 

def nfg0(e0):
    v,err=scipy.integrate.quad(cf0,0,e0)
    return v
dos=np.sqrt(En_mu)

f0=ff0(En_mu)
c_f0=cf0(En_mu)

f =  ff(En_mu)# Occupation number according to F-D statistics
df = np.gradient(f) # Take the first derivative of the F-D distribution.

c_f = cf(En_mu)  #np.sqrt(En_mu)/(np.exp( (factor*En_mu-1)/kt )+1) # Occupation number according to F-D statistics
c_df = np.gradient(c_f) # Take the first derivative of the F-D distribution.
n_fg=[nfg(x) for x in En_mu]
n_fg0=[nfg0(x) for x in En_mu]
n_f=[fg(x) for x in En_mu]

plot = figure(height=600, width=900, title="Change in the electron energy distribution with correction to chemical potential",
tools="pan,reset,save,wheel_zoom",
x_range=[0, 2], y_range=[0.0, 1.0])
plot.xaxis[0].axis_label='Energy (fraction of Ef)';
plot.yaxis[0].axis_label='dn/dE (Electron density per energy)';
r = plot.line(x=En_mu,y=f,legend_label="f(E) (Fermi-Dirac Distribution)", line_color = 'blue', line_width=2, line_alpha=0.6);
rr = plot.line(x=En_mu,y=dos,legend_label="g(E) (density of states)", line_color = 'brown', line_width=2, line_alpha=0.6);
rrr = plot.line(x=En_mu,y=c_f0,legend_label="g*f (the integrand)", line_color = 'red', line_width=2, line_alpha=0.6);
rrrr = plot.line(x=En_mu,y=c_f,legend_label="corrected g*f", line_color = 'green', line_width=2, line_alpha=0.6);
rrrrr = plot.line(x=[ef1,ef1],y=[0.0,2],legend_label="mu (chemical potential)", line_color = 'orange', line_width=2, line_alpha=0.6);
rrrrrr = plot.line(x=[ef,ef],y=[0.0,2],legend_label="Ef (Fermi energy)", line_color = 'black', line_width=2, line_alpha=0.6);
plot.legend.click_policy="hide"
show(row(plot),notebook_handle=True);


#Interactive part, letting user use the slider
def interactive_chem_consist(kt):
    '''
    :param k: wave number
    '''
    ef = 1
    kBT = 0.001 # Boltzmann constant measured interms of chemical potential
    E_min = 0.0001
    E_max = 3.0*ef
    factor = 1.0
    N = 100
    n0=50
    N1=20
    # Set up the variable space for 1D plot
    En_mu = np.linspace(E_min, E_max, N) #Energy measured in units of chemical potential
    En_l=np.linspace(E_min,ef,n0)
    En_r=np.linspace(ef,E_max,n0)
    
    
   
    ef1=ef
    nocc=2.0/3.0*(ef**(3.0/2.0))
    
    def nfg(e0):
        v,err=scipy.integrate.quad(cf,0,e0)
        return v
    def fg(e0):
        f1=lambda x: np.sqrt(x)
        v,err=scipy.integrate.quad(f1,0,e0)
        return v
    
    num_s=np.linspace(1,N1,N1)
    ef1_s=np.linspace(ef,ef,N1)
    for num in range(1,N1):
        cf=lambda x:np.sqrt(x)/(np.exp( (factor*x-ef1)/kt )+1) 
        c_f = cf(En_mu)  #np.sqrt(En_mu)/(np.exp( (factor*En_mu-1)/kt )+1) # Occupation number according to F-D statistics
        c_df = np.gradient(c_f) # Take the first derivative of the F-D distribution.


        #print(nfg(0.500))
        n_e=(1.5*nfg(E_max))**(2.0/3.0)
        
        ef1_s[num]=ef1
        ef1=ef1-0.2*(n_e-ef)
        
    ff0=lambda x:1.0/(np.exp( (factor*x-ef)/kt )+1)
    cf0=lambda x:np.sqrt(x)/(np.exp( (factor*x-ef)/kt )+1) 
    ff=lambda x:1.0/(np.exp( (factor*x-ef1)/kt )+1)
    cf=lambda x:np.sqrt(x)/(np.exp( (factor*x-ef1)/kt )+1) 
    
    def nfg0(e0):
        v,err=scipy.integrate.quad(cf0,0,e0)
        return v
    
    dos=np.sqrt(En_mu)
    
    f0=ff0(En_mu)
    c_f0=cf0(En_mu)
    
    f =  ff(En_mu)# Occupation number according to F-D statistics
    df = np.gradient(f) # Take the first derivative of the F-D distribution.

    
    c_f = cf(En_mu)  #np.sqrt(En_mu)/(np.exp( (factor*En_mu-1)/kt )+1) # Occupation number according to F-D statistics
    c_df = np.gradient(c_f) # Take the first derivative of the F-D distribution.
    n_fg=[nfg(x) for x in En_mu]
    n_fg0=[nfg0(x) for x in En_mu]
    n_f=[fg(x) for x in En_mu]

    r.data_source.data['y'] = f
    rr.data_source.data['y'] = dos
    rrr.data_source.data['y'] = c_f0
    rrrr.data_source.data['y'] = c_f
    rrrrr.data_source.data['x'] = [ef1,ef1]
    push_notebook();

    print("Electron Density, calculated at T = 0K (Ef taken as 1):",nocc)
    print ("Electron Density at finite temperature without fixing mu (chemical potential)",nfg0(E_max))
    print("Electron Density at finite temperature with fixed mu:",nfg(E_max))
    print("-------------")
    print("Corrected mu (chemical potential) as fraction of Fermi Energy=",ef1)
    print("Change in chemical potential d(mu):",np.abs((ef1-ef)/ef))
    print("-------------")
    print("Now to check with theory that change in chemical potential d(mu) = pi^2/12 * (kT/Ef)^2")
    print("12/pi^2 * d(mu)/(kT/Ef)^2 = ",np.abs(12/m.pi**2*(ef1-ef)/ef)/(kt/ef)**2, "[if second order approximation is accurate, should be 1]")

Note: Plot is interactive. If you feel like its too cluttered, click on legend to hide/show respective graph¶

In [3]:
style = {'description_width': 'initial'}
interact(interactive_chem_consist,kt= ipw.FloatSlider(min=0.01,max=0.2,step=0.005,description = 'kT as fraction of Fermi Energy', value=0.1, style=style, layout=Layout(width='500px')));

interactive(children=(FloatSlider(value=0.1, description='kT as fraction of Fermi Energy', layout=Layout(width…

Question:¶

So far we have finished the iterative calculation of chemical potantial variation. However, it is worth looking into the change of $E_f$ in reality.

Following section 2.2.3 in lecture notes, the chemical potential can be written (up to first term in Taylor expansion) as $$\mu=E_f \cdot [1-\frac{\pi ^2}{12}(\frac{k_B T}{E_f})^2]$$ Chemical potential changes from $E_f$ approximately as $(\frac{k_B T}{E_f})^2$. If the fermi energy in metals is in the scale of several eVs, would the correction be noticable at room temperature? You could use either the interactive slider or the approximation formula to answer.

In [ ]:
In [ ]: